LA-UR-07-7851 



Origin of the Cosmic Network: Nature vs Nurture 
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The large-scale structure of the Universe, as traced by the distribution of galaxies, is now being re- 
vealed by large-volume cosmological surveys. The structure is characterized by galaxies distributed 
along filaments, the filaments connecting in turn to form a percolating network. Our objective 
here is to quantitatively specify the underlying mechanisms that drive the formation of the cos- 
mic network: By combining percolation-based analyses with N-body simulations of gravitational 
structure formation, we elucidate how the network has its origin in the properties of the initial 
density field (nature) and how its contrast is then amplified by the nonlinear mapping induced by 
the gravitational instability (nurture). 

PACS numbers: 98.65.-r, 98.65.Dx, 98.80.-k, 98.80.Bp 



INTRODUCTION 



Observations of the large-scale distribution of galax- 
ies have been underway for several decades Going 
well beyond early "slice" views, recent observations re- 
veal a nontrivial three-dimensional structure 0, Q ■ Two 
key features of the structure are immediately apparent: 
(i) a considerable departure from local isotropy towards 
filament-like concentrations of galaxies and (ii) a ten- 
dency of the filaments to connect into a single perco- 
lating network spanning the entire region of observation. 
(A nice visualization of large scale structure as observed 
by the Sloan Digital Sky Survey can be found in Ref. 
A very similar picture also emerges from cosmological N- 
body simulations of the gravitational instability, the root 
cause of structure formation in the cosmological standard 
model (Fig. [J). 

Many statistics have been suggested and used to char- 
acterize cosmological structures. The galaxy two-point 
correlation function and its Fourier analog, the power 
spectrum, have been used extensively since the late 
1960s d, [g. As a way of studying the amplitude of den- 
sity fluctuations, these two-point statistics are observa- 
tionally robust and relatively straightforward to measure 
(in principle), as well as to predict theoretically (see, e.g., 
Ref. 0, and citations therein). However, being insensi- 
tive to phases, they do not contain shape information 
and cannot be used to probe the overall geometry and 
topology of the large-scale structure. The full hierarchy 
of n-point correlation functions, or the corresponding n- 
spectra, does carry complete information about the spa- 
tial distribution of galaxies. But n-point functions are 
difficult to measure as well as to predict. Furthermore, 
the desired information can be spread in a highly non- 
trivial manner over the space of n-point correlation func- 
tions. For these reasons, the use of higher-point statistics 
has been mainly restricted to the perturbative regime 
((((5/9_r/p)^)^/^ ^5 1) corresponding to smoothed fields 




FIG. 1: Demonstration of filamentary structure in N-body 
simulations. The projected density field is shown for a ACDM 
cosmology, with a box size of 64 ft-^Mpc. 



with a relatively large smoothing scale, i? ^ 5 h^^Mpc 
(p is the mean density and 5pR, the density fluctuation 
smoothed on a scale R as defined in Eqn. ([2]) below). 
Because of the lack of a single obvious way to proceed 
(and depending on the particular application in mind) 
many different statistics have been suggested for ana- 
lyzing large-scale structure observations and simulations. 
These include counts in cells, two- and higher order cor- 
relation functions i, the Euler characterrstic or equiva- 
lently the genus curve Q, global and partial Minkowski 
functionals the void probability function mini- 



mal spanning tree ll|, and many others 



In this work, we use percolation statistics to directly 
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address the issue of connectedness of superclusters of 
galaxies resulting in the emergence of the cosmic web. 
Our approach is attractive for two reasons: (i) The fact 
that simple scaling relations can be used to describe the 
percolation properties of cosmic density fields and (ii) 
that - as we show below - particle methods can also be 
conveniently folded into the analysis. 

The existence of the cosmic web in cosmological N- 
body simulations is certainly manifest (Fig. [T]). In search 
of a more analytic understanding of the dynamical ori- 
gin of the network topology, most theoretical appro aches 
start by invoking an idea due to Zel'dovich [l^, |l3| : ap- 
proximate the full evolution by a mapping from initial La- 
grangian (q) coordinates to the final Eulerian (x) state, 
via 



x(q,i) =q + i?(i)sflXq)- 



(1) 



Here x are the comoving coordinates, related to physical 
coordinates r as x = r/a{t), where a{t) is the scale factor 
describing the uniform expansion of the universe. The 
vector field Sij(q) is determined by the density fluctua- 
tions smoothed on a scale R at the initial time L , 



^fl(q, = [pR - p)/p = -D{ti)Vsi 



(2) 



and the monotonically growing function D(t) describes 
the growth of density fluctuations in the linear regime. 
Additional information present in the initial conditions 
may also be exploited, e.g., by employing the deforma- 
tion tensor djfc — —dsi/dqk and higher-order tensorial 
derivatives [IJ]. 

The dynamical mapping ([T]) is single- valued during the 
early (linear) phase of the evolution, becoming multi- 
valued after orbit-crossing (nonlinear phase). The ini- 
tial cosmic density field is a realization of a Gaussian 
random process, specified completely by the associated 
power spectrum, P{k). The choice of cosmological model 
fixes a{t), D{t), and P{k). We choose to set the smooth- 
ing of the initial perturbations to the present scale of 
nonlinearity, R « Rni, such that the linear rms density 
fluctuation a{Rni) = = 1 at the present epoch. 

This choice of the smoothing scale is based on the obser- 
vation (made in N-body simulations) that the large-scale 
density fleld is mainly determined by the linear power 
spectrum on the scales that have just become nonlin- 
ear, i.e., with k ^ R~l 15 1. It has also been shown that 



eliminating the perturbations on scales k ^ R~[ by trun- 
cation of the initial spectrum [l6l | or by smoothing with 
a Gaussian fllter [l3| removes structures at small scales 
but does not have much of an effect at large scales. 

One line of attack is to consider the evolution in a 
given realization of the initial conditions as a determin- 
istic map, and derive the properties of the network by 
studying generic singularities formed during nonlinear 
evolution 3i ld| using the deformation tensor, as well as 
higher order tensorial derivatives. Alternatively, one can 



employ a probabilistic approach, investigating network 
properties as determined by conditional multi-point cor- 
relation functions between certain parameters (e.g., the 
shear tensor) related to the deformation tensor at lin- 



ear density peaks [18|. These approaches have led to 



valuable insights, yet neither address the global proper- 
ties of the cosmic network directly. Both focus primarily 
on the structural building blocks (halos, filaments, pan- 
cakes) and may be loosely characterized as "local" , even 
though the characteristic scales of interest may reach tens 
of Mpc. 

In this paper, we take up the idea of identifying the 
cosmic network with percola ting regions in the density 
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field at the present time 



20( 1 . and apply percolation 



statistics as a quantitative global measure of the network 
structure 21| at both early and late times. In cosmo- 
logical applications of percolation, one studies the prop- 
erties of overdense and underdense excursion sets (ES), 
i.e., regions where the density is greater than some value 
{p > pc) or less {p < pc). Of particular interest are 
the connectivity properties of the excursion sets as the 
threshold pc is varied. 

Previous studies of percolation regions in the initial 
Gaussian and nonlinear fields ~ for two spatial dimen- 
sions - assumed the dynamics described by the Zel'dovich 
approximation and showed that the percolation region 
identified at the initial stage in the linear field became a 
part of the percolation region identified at the nonlinear 
stage [22|. The major reason for this conclusion was the 
continuity of the mapping described by the Zel'dovich 
approximation. Here, we investigate the relationship be- 
tween percolating regions in the initial (linear, Gaussian) 
and late-time evolved (nonlinear, non-Gaussian) density 
fields using N-body simulations of a ACDM cosmology. 
Our results establish the close - although not perfect - 
connection between the percolating regions in the linear 
and nonlinear stages of the evolution, as well as provide 
an understanding of the network structure based on per- 
colation concepts. 

Studies of the percolation properties of density fields 
obtained in cosmological N-body simulations with scale- 
free initial power spectra {P{k) oc /c") both in lin- 
ear and nonlinear regimes have been conducted previ- 
ously 23|, |2J] . In Ref . 23| a suite of N-body simulations 
in the Einstein-de Sitter universe (with n = —2, —1, 0, 1) 
was used to measure the percolation transition in the 
mass density field at four stages of dynamical evolution. 
It was shown that percolation transitions in overdense 
and underdense excursion sets experience shifts in op- 
posite directions from that of the initial Gaussian state 
(where they are obviously identical due to the symmetry 
of Gaussian fields with respect to the sign of the field) . In 
general, the percolation transition in the overdense phase 
occurs at lower than the Gaussian values while in the un- 
derdense phase it takes place at higher values of the vol- 
ume fraction (by volume fraction we mean the fraction 
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of the total volume in the corresponding phase). The 
dynamical evolution of the density field results in shrink- 
ing the percolating region in the overdense phase. At the 
same time the percolating region in the underdense phase 
gains considerably in volume. It is important to keep in 
mind that these are two independent statements since 
the overdense and underdense percolating regions have 
different boundaries in three-dimensional space. A sta- 
tistically significant trend was found relating the spectral 
index n to percolation: the greater the spectral index, the 
greater the delay in percolation as measured by larger val- 
ues of /p = (14s/14ot)at percolation, Or a Correspondingly 
lower density threshold for percolation, Sp/a. These re- 
sults were in qualitative agreement with Ref. [i^] which 
used a set of larger N-body simulations {n ~ —2, —1,0). 
In the latter work, the authors investigated the topol- 
ogy of the overdense and underdense excursion sets at 
the corresponding percolation thresholds in the density 
fields smoothed at various scales. They also measured 
percolation thresholds in the initial Gaussian fields and 
found the thresholds to be different for density fields with 
different power spectra. 



Our work differs from previous studies in several re- 
spects. First, we focus on the observationally realistic 
case of the ACDM model where the linear power spec- 
trum cannot be described by a single power law (previous 
studies of percolation in specific cosmological models do 
exist, e.g., Ref. however those models have since 

been ruled out observationally). Second, our simulation 
volume (both in terms of size and number of realizations) 
is significantly larger than considered previously. This 
helps to control finite size effects and reduces the statis- 
tical errors significantly. Finally, and most importantly, 
we consider a completely new dynamical problem - how 
does the (initial) linear percolating cluster map to the 
(final) nonlinear percolating cluster? By quantifying the 
accuracy of the mapping, we are able to distinguish be- 
tween two key aspects that underlie the structure of the 
cosmic web - an already present conspicuousness of the 
percolating region in the initial Gaussian field, which we 
refer to as nature, and the enhancement due to the grav- 
itational instability, which we term nurture. 



The rest of the paper is organized as follows: we first 
describe the cosmological models used in the study and 
then present results from percolation analysis of linear 
(Gaussian) density fields. After this we discuss the re- 
sults of the percolation analysis of the evolved nonlin- 
ear (non-Guassian) density fields (filtered with the same 
window as the linear density fields). In the last step, we 
study the accuracy of the mapping of the linear perco- 
lating cluster into the nonlinear percolating cluster, and 
end with a summary of the main results. 



PERCOLATION TRANSITION 
Cosmological Model 

The ACDM cosmological model adopted here is speci- 
fied by the following parameters: the dimensionless Hub- 
ble constant, h — i/o/(100 km/s/Mpc) ~ 0.7, dimension- 
less dark matter and baryon densities ^dm =0.259 and 
fib = 0.02 h^^, primordial spectral index n = 1, and a 
P{k) normalization set by as = (j{Rth = 8/i^^Mpc) = 
0.84, where Rth is the radius of the top-hat filter - the 
conventional measure of the amplitude of the initial den- 
sity fluctuations in cosmology. Our purpose is not to 
consider a model that is precisely the one given by cur- 
rent observations (which in any case would be a moving 
target), but is nonetheless close enough to have the same 
generic properties. 

Within this cosmology, we used two ensembles differ- 
ing in box size and filtering scale: 15 realizations gen- 
erated in a 341.3 /i~^Mpc box in linear dimension, fil- 
tered at i? = 1 /i^^Mpc (ACDMi) and 10 rcahzations 
in a 3413 /i^^Mpc box, filtered at i? = 10 /i^^Mpc 
(ACDM2). The study of the nonlinear dark matter den- 
sity field in ACDMi was carried out by using four statis- 
tically independent realizations of the initial conditions. 
(A nonlinear field smoothed with i? = 10 ft.~^Mpc ap- 
proximately corresponds to the linear field.) The nonlin- 
ear evolution is followed using an N-body particle-mesh 
code 26[, in a simulation cube of side 341.3 ft,~^Mpc. 



The number of particles equals the number of grid points 
Np = Ng ^ 512"^, in order to smoothly sample the den- 
sity field. (The force resolution is adequate given that 
here we are not interested in the small-scale distribution 
of matter.) In addition, we also considered two reference 
power-law models (n = — 2 and n = 4), for each of which 
10 realizations were generated. For all four ensembles, 
the filtering scale was 1.5 times the linear size of a spa- 
tial grid cell. 

Method 

Because the initial density field so strongly controls 
properties of the structure at late times 15, H, 17 1, we 
first discuss percolation in Gaussian fields with different 
P(fc), corresponding to four ensemble choices: two with 
ACDM P{k) (parameters as above) and two with power- 
law spectra P{k) oc fc" with n — —2 and 4 (all on a 512'^ 
grid). The power spectrum in the ACDM model is not a 
simple power law. Instead, in the limit of small fc, P{k) oc 
k while in the opposite limit, P oc fc^'^ln(fc/fcc), and at 
the present scale of nonlinearity (roughly), P oc k~^'^ 
(Fig. [2]). In the not so distant past the scale of nonlin- 
earity was smaller and the effective slope was steeper, ap- 
proximately n = — 2; in principle this evolutionary stage 
is observable. The choice of n = 4 corresponds to the 



4 



100000 



10000 



1000 




100 



0.1 
k [h/Mpc] 

FIG. 2: Power spectra of two ACDM models with smooth- 
ing scales, R — 1 and 10 h~^Mpc, and box-sizes 343.1 and 
3431 h'^Mpc, respectively. The unsmoothed ACDM power 
spectrum and P{k) oc are also shown for reference. 



SO called 'minimal power spectrum' on large scales [27 1 
(this power spectrum would arise from initial conditions 
which have power only on small spatial scales). 

Percolating sites are identified using the "friends-of- 
friends" algorithm on a cubic mesh (only the eight clos- 
est neighbors were considered as immediate friends) for 
the (thresholded) density field filtered on a scale R with 
a Gaussian window. As the density threshold is reduced 
from a high starting value, the volume in the overdense 
excursion set increases monotonically until percolation 
occurs, therefore the volume fraction of the excursion set 
itself can be used as a proxy for the density threshold. 
(The same argument also holds for the underdense excur- 
sion set.) For the linear density field, percolation curves 
were averaged over overdense and underdense excursion 
sets exploiting the exact statistical symmetry of Gaus- 
sian fields. The power spectra in the unsmoothed ACDM 
model are shown in Fig. [2] along with the power spectra 
in all realizations of the ACDMi and ACDM2 models; a 
power-law spectrum with P{k) oc is also shown for 
reference. 



Linear Stage 

We first focus on the initial stage of structure forma- 
tion, i.e., any epoch after decoupling of baryonic matter 
and radiation with only linear density fluctuations on the 
scales of interest. For the chosen resolution scales (set by 
the particle number and box size), an initial redshift of 
Zi = 50 comfortably satisfies this criterion. 

In order to carry out our percolation analysis we used 
the volume fraction of the excursion set rather than 
the density itself as the relevant parameter. The re- 
sulting data - the percolation curves - turn out to be 
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FIG. 3: Percolation transitions in Gaussian density fields. 
Results for power-law spectra P{k) oc fc" with n = —2 (ex- 
treme left) and n = 4 (extreme right), ACDM - the two inte- 
rior curves - smoothed at 7? = 1 (left) and i? = 10 /i"^Mpc 
(right). Smooth curves are fits to the percolation ansatz ((3|; 
error bars are la (statistical). 



TABLE I: Values for the percolation ansatz ([Sjl parameters {A 
and v) for Gaussian fields, the percolation density threshold 
Sp corresponds to /p, the volume fraction at percolation. 



Model 


/p 


A 


V 


5j,/a 


n=4 


0.157 


0.61 


0.38 


1.006 


ACDM2 


0.111 


0.66 


0.51 


1.22 


ACDMi 


0.089 


0.75 


0.62 


1.36 


n=-2 


0.072 


0.89 


0.76 


1.46 



well-described by the single-variable percolation scaling 
ansatz 1281 



/i = ^(/es - /p)" at /es>/p, 



(3) 



where /es is the volume fraction, i.e., the fraction of the 
volume in the excursion set (/es = VEs/Kot), A, the 
volume fraction of the largest region, and, /p, the value 
of the volume fraction when the percolation transition 
occurs. In Fig. [3l the quantity, /i/Zes, as numerically 
determined, is plotted as a function of /es, and undergoes 
a sudden growth from zero to unity at the percolation 
transition. As is evident, the fitting form ^ provides a 
good fit to the data for all four models. The parameters 
of the fits are given in Table [D 

The data shown in Fig. [3] clearly demonstrate that per- 
colation occurs at lower volume fractions with growth 
of power on large scales, in qualitative agreement with 
Ref. [iJ]. The corresponding density threshold reaches 
its (lower) hmiting value {5p = cr,fn = 0.157), fol- 
lowing a conjecture due to Ziman [29| who suggested 
that the percolation transition should occur exactly at 
the 6 — a level, for Gaussian fields with dominating 
small-scale power. At the other extreme considered here. 
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FIG. 4; The function A^(fc) = k^P{k) plotted for the power 
spectra shown in Fig. (2] The power law spectrum with n = —2 
is smoothed with the same window as the ACDMi model. The 
normalization of A'^(fc) is not important since it determines 
only the variance of the field. 



P{k) oc fc^^, where the field is dominated by large-scale 
power, /p = 0.072, and corresponds to the highest value 
of the density threshold 5p/a = 1.46. This is slightly 
lower than 5p/a = 1.6±0.1 as found in Ref. ^24j]. The 
minor disagreement may be due to how the percolation 
threshold was determined in the two cases, as well as 
due to finite- volume limitations. [It is also worth noting 
that our set of fields includes models with nonpower law 
P{k)] Each parameter in Table [T] varies monotonically 
with increase of power on large scales: /p decreases while 
all the rest increase. 

From the results of Fig. [3] it is clear that the perco- 
lation transition becomes more gradual as the (effective) 
spectral index n decreases, deviating from the percolation 
scaling ansatz at small /i / /es ■ Therefore the fit given by 
Eqn. ([3]) becomes less accurate and /p determined by this 
method becomes less reliable. This suggests that even 
larger volumes than the ones used here may be necessary 
to study such cases. Other factors that may affect the ac- 
curacy of the determination of the percolation threshold 
are discreteness, accuracy of computing the volume of the 
excursion set, filtering, and the properties of the largest 
cluster. As the present study is not devoted to an ac- 
curate calculation of percolation thresholds in Gaussian 
fields, we leave these questions to future work. 

For scale- free power spectra, P{k) oc fc", the depen- 
dence of the percolation transition in Gaussian fields on 
the power spectrum was observed and qualitatively dis- 
cussed in Ref. [2^ and then quantified in Ref. 24 1. Here, 
along with two types of scale-free fields, we also study 
the Gaussian field corresponding to the ACDM model 
smoothed at two different scales (Fig. [2]). As mentioned 
earlier the large-scale structure at the nonlinear stage 
is mainly determined by the initial power spectrum on 



scales with k ^ R and therefore can be ap- 

proximately characterized by the effective slope at this 
scale. It is thus interesting to investigate how well the 
effective slope at the scales contributing most to the vari- 
ance of the field determines the percolation transition in 
the linear regime. 

As a first candidate, we study A'^{k) = k^P{k), which 
represents the contribution of power to the variance of 
the density field, per unit logarithmic interval in k. In 
particular, the position of the maximum of A^(A:) corre- 
sponds to the wavelength with maximum contribution to 
the variance of the field. Figure H shows that the maxi- 
mum contribution arises from k « 0.65 and 0.09 /iMpc~^ 
for the ACDMi and ACDM2 models respectively. The 
initial i.e., unsmoothed spectrum has the corresponding 
effective slopes ricg = dlogP/dlogk « —2.2 and —1.4 at 
these wavenumbers. 

Referring to Table HI two subsets of the models rep- 
resented there - one subset consisting of the power law 
model with n = — 2 , ACDM2 with rics = —1.4, and 
the power law model with n = 4, the other consisting of 
ACDMi with Ucs = -2.2, ACDM2 with rics = -1-4, 
and the power law model with n = 4, both demon- 
strate monotonic growth of the percolation volume frac- 
tion with increase of effective slope: fp — 0.07, 0.11, 0.16 
and fp — 0.09, 0.11, 016 respectively. However, perco- 
lation in the ACDMi model with Ues — —2.2 formally 
takes place at a higher value of the volume fraction, fp 
— 0.09, than for the power law model with n — —2, 
where fp = 0.07. Regardless of the limited accuracy in- 
herent in measuring the values of fp from the percolation 
ansatz ([3]), it does appear rehable to conclude directly 
from the data points of Fig. [3] that the percolation tran- 
sition in the power law model with n = —2 occurs at a 
lower volume fraction then in the ACDMi model. Addi- 
tionally, the data in Fig.[2]show that the effective slope of 
the unsmoothed ACDM model is slightly more negative 
than n = —2 at fc « 0.9 hMpc~^ (where A^(/c) peaks for 
the n — —2 power spectrum) . 

An alternative possibility is to use the power per equal 
interval in k rather than in log/c, i.e., E{k) = k^P(k) in- 
stead of A^(fc) — k^P{k); the former choice being widely 
used in the theory of turbulence [sO] ■ Again, by finding 
the maximum of E{k) and then the effective slope in the 
unsmoothed spectra at these points, one finds that all 
four models listed in ascending order of fp (i.e., n = —2, 
ACDMi, ACDMi, and n = 4) correspond to monoton- 
ically increasing effective slopes, — -2, -1.97, -1.1, 
and 4. The difference in rieff between the power law 
model with n = —2 and the ACDMi model is obviously 
marginal, but the effective slopes at the maxima of E{k) 
do appear to correspond better to the volume fractions at 
percolation. The maxima of E{k) are obviously shifted to 
smaller k with respect to the maxima of A^(fc), therefore 
E{k) takes into account power on slightly smaller scales 
than the scales of the maxima of A^(fc). A more detailed 
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evaluation of this issue is a topic for further study and 
beyond the scope of the current work. 

Despite the quahfications in the above discussion, the 
ACDM percolation curves do display trends qualitatively 
similar to those seen in the bracketing power-law cases. 
Filtering the density field at smaller scales corresponds 
to a greater negative effective spectral slope than that 
due to filtering at larger scales (Fig. [2]) , the effect of the 
moving cutoff in /c-space leading therefore to a higher 
density threshold and lower volume fraction at percola- 
tion, as demonstrated by the two interior curves in Fig. [31 
In particular, this means that the structure formed ear- 
lier at greater redshifts was formed from the Gaussian 
field percolating better than the field corresponding to 
the current structure. 

The dependence of the percolation transition param- 
eters on the power spectrum of a Gaussian field shows 
that such fields can differ not only locally (e.g., in the 
shape of peaks, etc. [11]) but also globally. The higher 
percolation threshold, Sp, or equivalently lower value of 
/p, indicates that Gaussian fields with relatively more 
power on large scales have a greater degree of connected- 
ness. Therefore, the relatively large negative slope of the 
linear power spectrum at the scale of nonlinearity in the 
ACDM cosmology is an important factor determining the 
origin of the cosmic web. Because this feature is already 
present in the initial conditions, later to be amplified by 
the gravitational instability, we refer to it as the 'nature' 
factor influencing the formation of the cosmic web. 
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FIG. 5: ACDMi percolation transitions (filter scale set to 
R = 1 h^^Mpc). Nonlinear stage: overdensity (blue/left) 
and underdensity (red/right) excursion sets at z = 0; linear 
(Gaussian) stage (green/middle) nt z — 50. Solid curves are 
percolation ansatz fits (cf. Eqn. |3]and Table HI)) . 



TABLE II: Values for the percolation ansatz (Eqn. [Sj pa- 
rameters for underdense and overdense excursion sets at the 
nonlinear stage of evolution, the percolation density threshold 
Sp corresponds to /p. 



Model 


/p 


A 


V 


Sp 


NL underdense 


0.228 


1.80 


0.76 


-0.80 


NL overdense 


0.035 


0.73 


0.75 


3.31 



Nonlinear Stage 

Next we turn to a percolation transition analysis of 
the nonlinear, and therefore non-Gaussian, density field 
at z = for the ACDMi model. As shown in Fig. [5l 
the percolation threshold for the overdense excursion set 
at the nonlinear stage (as measured by the volume frac- 
tion) is remarkably lower than that for the initial Gaus- 
sian field: /p(z = 0) = 0.035 while /p(z = 50) = 0.089. 
This result is as expected from the Zel'dovich approxi- 
mation [i^]. The evolution of the density field results 
roughly in compression of overdense regions and expan- 
sion of underdense regions. This qualitatively correct 
feature is complicated by the fragmentation of the dense 
regions into gravitationally bound halos that might make 
percolation more difficult. However, the relative strength 
of both effects strongly depends on the power spectrum 
of the initial Gaussian field. This study suggests that de- 
spite the strong fragmentation of the linear percolating 
region, the overall spatial distribution of the fragments 
markedly preserves the topology of the linear percolating 
region, at least in the important case of ACDM. 

Compared to the linear case, the density threshold for 
percolation in the nonlinear regime is significantly larger. 
(The mean values and standard deviations are obtained 



from four realizations of the nonlinear stage and thirty re- 
alizations of the linear Gaussian field.) A visual compar- 
ison can be made by comparing the nonlinear and linear 
percolating clusters. These are shown in Figs. [H] (the blue 
region is the linear percolating cluster) and [7] (the yellow 
region is the nonlinear percolating cluster). At the same 
time the total mass in the nonlinear percolating cluster is 
about 15 times greater than that in the corresponding lin- 
ear cluster: the mass fraction Mi « 0.25 at z = 0, while 
Ml « 0.018 at z = 50. The rightmost curve in Fig. [S] 
shows that the percolation transition in the underdense 
excursion set occurs at considerably greater volume frac- 
tion than in the initial Gaussian field: /p(z = 0) = 0.23 
while /p(z = 50) = 0.089. 

The parameters of the fitting curves are given in Ta- 
ble |TT1 It is worth noting that the percolation transitions 
in all four tested Gaussian fields are significantly greater 
than the percolation threshold - measured by volume 
fraction - for the over-dense phase but smaller than that 
for the under-dense phase of the nonlinear density field 
generated by the gravitational instability. Expanding the 
study of percolation to n — > —3 is relatively difficult, 
therefore whether the above result holds in this limit re- 
mains to be seen. It is also worth pointing out that the 
values of the exponent v are similar for both overdense 
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and underdense excursion sets, approximately coinciding 
with the value of v obtained for the scale-free Gaussian 
field with n = — 2. The amplitude A is different in all 
three cases, however. 

MAPPING ACCURACY 

Standard Eulerian perturbation theory at linear order 
does not move particles. However, if one takes the perco- 
lating cluster in the linear density field, smooths it at the 
scale of nonlinearity, and maps it using the Zel'dovich ap- 
proximation then in effect one has a (linear theory) pre- 
diction for the nonlinear percolating cluster. Intuitively, 
one may expect that the percolating cluster in the linear 
field should approximately map into the exact nonlinear 
percolating cluster, but the accuracy of this proposition 
has not been directly tested. In this section, we focus on 
whether the linear percolating cluster indeed maps onto 
the nonlinear percolating cluster. 

We begin by considering the mapping from the initial, 
linear stage of the evolution [z = 50), to the final non- 
linear stage (z = 0) , the forward map. This requires two 
additional steps. First, we choose a fiducial value for 
the density threshold, S''^\ to carry out the percolation 
analysis. We set this to correspond to the point where 
the largest percolating region occupies 50% of the excur- 
sion set volume, i.e., /{'^V/es = 0.5. With this choice, 
the percolation region is well developed, yet not com- 
pletely dominating the excursion set. This definition is 
different from that used in the previous analysis (which 
followed the conventional practice in percolation theory) . 
The reason for this change is dictated by the fact that at 
the (conventional) percolation transition, the percolating 
cluster occupies an infinitesimal volume (/i//es 0). 
This means that in a finite box its volume is much smaller 
than the volume of the box and is thus subject to very 
large fiuctuations. Our choice makes the definition of 
the percolating cluster more robust and less subject to 
the precise "turn-on" of percolation (Cf. Figs. [3] and O. 

The second step involves multiple back and forth 
switching from density fields defined on an invariant spa- 
tial grid to particles subject to dynamic motion. Al- 
though the percolation analysis can be carried out on a 
set of points it is numerically more efficient on a reg- 
ular grid. In this study we used a very fast algorithm 
that detects the percolation transition by performing the 
cluster analysis on the grid [2^. In particular, it defines 
the percolating cluster as a set of grid sites. In order 
to study the forward mapping of the initial percolating 
cluster (identified as a set of grid sites), we first need to 
determine the set of particles that are associated with it. 
Following the motion of these selected particles will allow 
us to study how the initial percolating region maps into 
the nonlinear stage. 

As a definition, the initial percolating set is taken to 



be all those particles lying within a distance smaller than 
half a grid unit (along all three orthogonal axes) from 
any grid site that belongs to the initial percolating re- 
gion. After mapping particles to the final nonlinear state 
we find all the sites on the grid that satisfy the very same 
criterion. The set of these new sites at z = approxi- 
mates the map of the initial percolating cluster and can 
be analyzed by the same fast method. We stress that 
this forward mapping of the linear percolating cluster 
comprises only a relatively small fraction of particles in 
the nonlinear percolating cluster because a fundamental 
feature of the nonlinear evolution is the crossing of orbits 
and the resulting formation of multi-stream flow regions. 
The forward mapping of the percolating cluster found 
in the linear density field forms only a fraction of those 
streams, and many others come from regions beyond it. 
In a similar manner, by starting with the nonlinear perco- 
lating cluster at z = 0, we introduce a backward mapping 
that finds all sites at z = 50 that were mapped to z = 
and formed the nonlinear percolating cluster. 

With the particle-based definition of percolation ap- 
plied to the ACDMi cosmology, we find that the particles 
in the initial percolating set make up 4.6% of the total 
number of particles in the excursion set - quite close to 
the volume fraction of the percolating region as defined 
on the grid (0.045). The small excess is due to the higher 
than mean density in the percolating region: the two def- 
initions of a percolating region (on grid and by particles) 
are not identical. But we do expect them to be very sim- 
ilar during the linear stage of the evolution. To verify 
this expectation we go back from the selected particles 
to the grid by identifying every grid site closer than half 
a grid unit to any particle from the percolating set. We 
then perform a cluster analysis on these sites. Ideally we 
would find exactly the same grid sites as in the initial 
percolating cluster. The actual results show good but 
not perfect correspondence between the two procedures. 

Figure [6] illustrates the forward mapping of the initial 
percolating set (blue) into the nonlinear stage (yellow). 
Both structures are displayed in comoving coordinates 
to scale out the uniform Hubble expansion. To avoid 
excessive projections, every tenth particle is shown in a 
70 h~^Mpc thick slab cut out of the simulation cube. The 
mapping obviously results in a large reduction of the vol- 
ume as expected, therefore the nonlinear (yellow) cluster 
is much smaller by volume than the corresponding linear 
cluster (blue). A closer look also suggests that the ini- 
tially connected structure at z = 50 becomes fragmented 
after nonlinear mapping to z = 0. We discuss this in 
more detail below. 

In the first step of the analysis, we compute the fraction 
of particles from the initial percolating set that end up 
in the final percolating region, as defined on the grid at 
various thresholds, and shown for the fiducial S^^^ — 2.76 
in Fig. [7] in yellow. (We determine whether a particle be- 
longs to the percolating region at z = using the same 
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FIG. 6: Forward mapping of the linear percolating cluster. The mapping takes the percolating region identified in the linear 
regime at 2 = 50 (blue) to the final forward- mapped region at 2 = (yellow, see text). Note the fragmentation of the initial 
percolating region after the mapping. For visual clarity, a 70 h~^Mpc thick slab cut out of the full simulation cube is shown. 




FIG. 7: Backward mapping of the percolating region identified in the nonlinear regime at 2 = (yellow) to the initial stage 
at 2 = 50 (blue). In contrast to the case of the forward map, the backwards map takes a percolating region into another 
percolating region. The slab is the same as in Fig. [6l 
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criterion of proximity of particles to grid sites as above, 
filtering the nonlinear density field also at R=l h~^Mpc.) 
The fraction of particles from the initial percolating set 
transported into the final percolating region is plotted 
as a function of the ratio fi/fEsiz — 0), character- 
izing the size of the percolating cluster at z = 0, in 
Fig. [8l This fraction grows monotonically from about 
40% at /i//es(z = 0) « 0.22 to more than 95% at 
/i//es(-z ~ 0) ~ 0.8 reaching about 80% at the fidu- 
cial value /i //es(-z = 0) = 0.5. The point of this result 
is that the progenitor of the "cosmic web" defined as the 
percolating cluster of the initial Gaussian density field, 
is a fair, albeit not perfect, "backbone" of the web at a 
later, nonlinearly evolved stage. 

It is also of interest to investigate the percolation prop- 
erties of the initial set of percolating particles at the final 
stage of evolution. Using the criterion described above, 
we generated the set of grid sites at z = close to the se- 
lected particles (those representing the percolating clus- 
ter at z — 50) after they had moved to their final po- 
sitions at z = 0. We find that these sites comprise less 
than 0.003 of the total volume. This should be compared 
with their initial volume occupation fraction of 0.045 at 
z = 50, resulting in a compression factor of about 15 
on average. This compression reflects the fact that over- 
dense regions tend to collapse and therefore more than 
one particle becomes associated by the proximity cri- 
terion to the same site at the nonlinear stage. This is 
clearly seen in Fig. [6l By z = 0, the initially connected 
region fragments into large numbers (~ 2.5 x 10^) of iso- 
lated regions and none of them percolates. The largest 
such region has a volume of ~2000 h~^Mpc^. Although 
this volume is quite large by astronomical standards (al- 
beit much smaller than that set by the scale of homo- 
geneity which has a linear dimension of ^^70 /i~^Mpc), it 
comprises only 0.017 of the volume occupied by the parti- 
cles from the initial percolating set after they are mapped 
to z = 0. These numbers illustrate the prolific fragmen- 
tation of the linear percolating cluster in the course of 
the nonlinear mapping. 

We now take the percolating region at z = 0, shown 
in Fig. [7] in yellow, back to the linear stage - the back- 
ward map (blue region in Fig. [7] with every twentieth 
particle plotted). As an example, we select the perco- 
lating region built from the smoothed density field at 
df = 2.76 corresponding to fi^/hs{z = 0) = 0.5. This 
occupies 8.4 x 10^ h^'^Mpc' or 0.02 of all grid sites or 
equivalently about 0.02 of the volume of the box. The 
remaining half of the excursion set consists of ~ 2 x lO'' 
isolated regions. It is worth stressing that approximat- 
ing the percolating region at the nonlinear stage by the 
particle set as described above, we find that the two rep- 
resentations (particle vs. grid) differ considerably more 
than they do in the linear regime. This is because the 
particle representation has a significantly higher resolu- 
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FIG. 8: Fraction of the particles from the initial percolating 
set {z = 50) in the percolating region at 2; = as a function 
of /i//es at 2: = 0. 

tion than the smoothed density field on the grid. In any 
case, the backward map must be applied using the par- 
ticle representation; inverse mapping the particles from 
the largest percolating set back to the linear stage, we 
find that they occupy about a quarter of the total vol- 
ume and also form a percolating region (in contrast to 
the case of the forward map). 

SUMMARY 

We now summarize our key results. We began by con- 
firming - with some minor discrepancies - earlier findings 
of the dependence of the percolation threshold in Gaus- 
sian random fields on the power spectrum. However, the 
purpose of our study was not the precise measurement 
of the percolation thresholds for Gaussian fields but to 
establish the existence of strong percolation in the ini- 
tial/linear density field in the ACDM cosmological model 
- the nature factor in the formation of the cosmic web. 
This factor is determined exclusively by the density (or 
density contrast) field, 5 = {p — p)/p in Eulerian linear 
theory. 

We found that the nonlinear mapping from the linear 
stage of cosmic evolution (z = 50) to the final nonlinear 
stage at the current epoch (z = 0) results in a consid- 
erable amplification of percolation features: the percola- 
tion volume fraction reduces from /p(z = 50) « 0.09 to 
/p(z = 0) « 0.035, the half-fill volume of the percolating 
cluster at the fiducial percolation volume fraction reduces 
from f[^\z = 50) « 0.05 to f[^\z = 0) « 0.02, while 
its mass greatly increases from vrSi\z = 50) w 0.025 to 
rrv{ {z = 0) w 0.25. All these parameter values indi- 
cate that the cosmic web becomes considerably thinner 
and more conspicuous than its progenitor in the linear 
density field. We call this dynamical infiuence of the 
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gravitational instability, the nurture factor in the forma- 
tion of the cosmic web because it was not present in the 
initial field and results entirely from the mapping itself. 
Conventional linear theory does not consider mapping at 
all; the Zel'dovich approximation suggests that the initial 
displacement field, s = — 9$iin/(9q smoothed at the cur- 
rent scale of nonlinearity determines the dominant fea- 
tures of the mapping. 

We conclude that the robust contrast of network struc- 
ture - so characteristic of N-body simulations - is de- 
termined by a combination of two factors: the already- 
present conspicuousness of the percolating region in the 
initial Gaussian field (nature) and the effects of nonlinear 
mapping [nurture). The first is quantified by the rela- 
tively high density threshold [S'^^^ = 1.27a) and low vol- 
ume fraction of the excursion set (/^•^•' « 0.1) at the per- 
colation transition. This is directly related to the charac- 
ter of the linear power spectrum - a relatively large neg- 
ative effective slope in the relevant range of scales from 
k = 0.1 to 1 hMpc^^. The second factor is measured by 
the further decrease of the volume fraction (/'■•^■' ~ 0.05) 
at percolation. The percolation threshold in the density 
field filtered on i? = 1 h'^Mpc is S^f) = 2.76 with the 
percolating region containing about a quarter of the to- 
tal mass. We conclude that the nurture factor arising 
from the dynamical mapping itself is more important for 
explaining the conspicuousness of the cosmic web. 

The topology of the web can also be quantified by the 
Euler characteristic, but is not equivalent to percolation 
statistics [sij. In particular, the density of the Euler 
characteristic, %, as a function of the level, — S/as 
has the universal shape, n^^ — N^{v'^ — 1) exp(— :/^/2) 
for all Gaussian fields regardless of the power spec- 
trum ^33|. The power spectrum determines only its 
amplitude, = 2((fc2)/3]3/2/(27r)2, where (fc)^ = 

/ PPik)Pdk/ J P{k)Pdk 0]. Thus, the standard di- 
agnostics used in cosmology for characterization of the 
structure: "meatball shift", Ai/, "number of voids", Ay, 
and "number of clusters", Ac (see e.g., Ref. [3J]) de- 
signed for studies of non-Gaussian features of the struc- 
ture are unable to detect the 'nature' factor in structure 
formation. 

We note that the mapping of the linear percolating 
cluster into the nonlinear stage is not flawless. The per- 
colating region at the nonlinear stage taken at the half- 
fill percolation volume fraction contains about 80% of the 
percolating region found in the linear density field, also 
taken at the half-fill percolation volume fraction. This is 
close, though not perfect, correspondence. The half-fill 
linear percolating cluster makes up about 5%, while the 
nonlinear one about a quarter, of the total mass, i.e., five 
times greater. The initial percolating region fragments in 
the course of evolution into a large number (more than 
2.5 X lO'' in the current simulation) of isolated clumps by 
z = 0, none of which percolates at z = 0. 

A comparison of the percolation properties of the 



galaxy distributions in the ACDM model and galaxy red- 
shift surveys requires mock redshift catalogs obtained 
from simulations covering volumes similar to the 2dF or 
SDSS surveys. An additional theoretical complication 
is related to the fact that real galaxy distributions are 
available only in redshift space. Redshift space is fun- 
damentally anisotropic in the sense that the statistical 
properties of the large-scale structure are different along 
the line of sight and in the orthogonal directions [s^. 
This may require additional refining of the technique dis- 
cussed here for its analysis. However, some observational 
studies have already been carried out and have succeeded 
in measuring a significant signal of filamentarity in red- 
shift catalogs [36| . 
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